rm(list=ls())
library(foreign)
library(data.table)
library(RColorBrewer)
library(grid)
library(ggmap)
library(maptools)
library(dplyr)

# shapefiles 

cz <- readShapePoly("CZ.shp")
cz.df <- fortify(cz)
cz.df <- data.frame(cz.df, cz@data[cz.df$id, ])
cz.df <- cz.df[cz.df$long < -50 & cz.df$long > -135 & cz.df$lat < 50 & cz.df$lat > 22, ]

names(cz.df)[8] <- "czone"

perc.rank <- function(x) trunc(rank(x))/length(x)

# reading czone-level data

cz_taa <- read.dta("kimpelc_czone.dta")
cz_taa$mean_l_TAAdlrs_wrks_cert_pc_qt <-  perc.rank(cz_taa$mean_l_TAAdlrs_wrks_cert_pc)
cz_taa$mean_l_trans_taaimp_pc_qt <-  perc.rank(cz_taa$mean_l_trans_taaimp_pc)
cz_taa$difference <- cz_taa$mean_l_TAAdlrs_wrks_cert_pc - cz_taa$mean_l_trans_taaimp_pc
cz_taa$difference_break <- cut(cz_taa$difference, breaks = c(-20000, -2000, -1500, -1000, -500, 0, 500, 1000, 1500, 2000, 20000))
cz_taa $success_rate <- cz_taa $certified_petitions_n/cz_taa $petitions_n

cz.df<- merge(cz.df, cz_taa, by = intersect(names(cz.df), names(cz_taa)))

# map base 

mapPlot <- ggplot(cz.df) + theme_bw() + theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank() ,axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), legend.position="bottom",    legend.key.width = unit(1, "cm"), legend.text = element_text(size = 8))

# Figure 1-A TAA Spending Based on Unemployment Claims

pdf("Figure1-A.pdf",width=7,height=5)

adh_map <- mapPlot + geom_polygon(data=cz.df, aes(x = long, y = lat, group = group, fill = cz.df$mean_l_trans_taaimp_pc_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "Average TAA Spending per 1,000 Workers", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0), labels=c("6", "661", "1,106", "2,037", "11,846")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(adh_map)
dev.off()

# Figure 1-B TAA Spending Based on TAA Certified Workers

pdf("Figure1-B.pdf",width=7,height=5)

adh_map <- mapPlot + geom_polygon(data=cz.df, aes(x = long, y = lat, group = group, fill = cz.df$mean_l_TAAdlrs_wrks_cert_pc_qt), lwd = 1/10, colour="lavenderblush4") +  scale_fill_continuous(name = "Average TAA Spending per 1,000 Workers", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "50", "794", "2,244", "23,299")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(adh_map)
dev.off()

# reading czone-by-decade data

cz.response <- read.dta("kimpelc_czone_t2.dta")
cz.response <- subset(cz.response, year!= 2007)

cz.response$exposure_qt <- perc.rank(cz.response$l_tradeusch_pw)
cz.response$response_qt <- perc.rank(cz.response$response)
cz.response$taa_qt <- perc.rank(cz.response$petitioningwrks_pw_mean)
cz.response$cert_response_qt <- perc.rank(cz.response$response_cert)
cz.response$response_qt9010 <- perc.rank(cz.response$response9010)

cz.response.1990 <- subset(cz.response, t2 == 0)
cz.response.2000 <- subset(cz.response, t2 == 1)

cz.change <- merge(cz.response.1990, cz.response.2000, by = c("czone"))

cz.change$exposure_delta <- cz.change$l_tradeusch_pw.y - cz.change$l_tradeusch_pw.x
cz.change$exposure_delta_qt <- perc.rank(cz.change$exposure_delta)

cz.change$response_delta <- cz.change$response.y - cz.change$response.x
cz.change$response_delta_qt <- perc.rank(cz.change$response_delta)

cz.change$taa_delta <- cz.change$petitioningwrks_pw_mean.y - cz.change$petitioningwrks_pw_mean.x
cz.change$taa_delta_qt <- perc.rank(cz.change$taa_delta)

cz.change$response_delta9010 <- cz.change$response9010.y - cz.change$response9010.x
cz.change$response_delta_qt9010 <- perc.rank(cz.change$response_delta9010)

cz.change <- cz.change[,c("czone", "response_delta", "exposure_delta", "taa_delta", "response_delta_qt", "exposure_delta_qt", "taa_delta_qt", "response_delta9010", "response_delta_qt9010")]

cz.df.1990 <- merge(cz.df, cz.response.1990, by = intersect(names(cz.df), names(cz.response.1990)))
cz.df.2000 <- merge(cz.df, cz.response.2000, by = intersect(names(cz.df), names(cz.response.2000)))

cz.df.change <- merge(cz.df, cz.change, by = intersect(names(cz.df), names(cz.change)))

# Figure 2 TAA Petitioners by Import Exposure per Worker in the 1990s (A) and 2000s (B) -- 2000-2007

pdf("Figure2-A.pdf",width=7,height=5)

response_1990 <- mapPlot + geom_polygon(data=cz.df.1990, aes(x = long, y = lat, group = group, fill = cz.df.1990$response_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "TAA Response to Import Exposure", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.066", "0.285", "0.772", "770.04")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_1990)
dev.off()

pdf("Figure2-B.pdf",width=7,height=5)

response_2000 <- mapPlot + geom_polygon(data=cz.df.2000, aes(x = long, y = lat, group = group, fill = cz.df.2000$response_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "TAA Response to Import Exposure", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.066", "0.285", "0.772", "770.04")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_2000)
dev.off()